import warnings 
warnings.filterwarnings("ignore")

from IPython.core.display import display, HTML
from IPython.display import Image, display, HTML

HTML('''<script>
code_show=true; 
function code_toggle() {
 if (code_show){
 $('div.input').hide();
 } else {
 $('div.input').show();
 }
 code_show = !code_show
} 
$( document ).ready(code_toggle);
</script>
<form action="javascript:code_toggle()"><input type="submit" value="Click here to toggle on/off the raw code."></form>''')

# For Data Prep
import numpy as np
import pandas as pd
import ast
import sqlite3
from sqlalchemy import create_engine
from sklearn import preprocessing
from sklearn.preprocessing import StandardScaler
from tqdm import tqdm
from collections import Counter
import datetime
from sklearn.feature_extraction.text import TfidfVectorizer
from sklearn.preprocessing import MultiLabelBinarizer
from wordcloud import WordCloud, STOPWORDS
from PIL import Image
from tqdm import trange

# For data scaling
from sklearn.preprocessing import StandardScaler
from scipy.stats import norm

# For Dimensionality Reduction
from sklearn.decomposition import PCA, TruncatedSVD

# For Representative-Based Clustering
from sklearn.cluster import KMeans
from pyclustering.cluster.kmedoids import kmedoids
from sklearn.base import clone
from sklearn.metrics import calinski_harabasz_score, silhouette_score
from sklearn.metrics import davies_bouldin_score
from scipy.spatial.distance import euclidean, cityblock

# For Heirarchical
from scipy.cluster.hierarchy import dendrogram
from scipy.cluster.hierarchy import fcluster
from fastcluster import linkage

# For Density-Based Clustering
from sklearn.neighbors import NearestNeighbors
from sklearn.cluster import DBSCAN
from sklearn.cluster import OPTICS, cluster_optics_dbscan

# For Gaussian Mixture Clustering
from sklearn.mixture import GaussianMixture as GMM

# For Recommender System
from surprise import (Reader, Dataset, KNNWithMeans, SVD)
import recmetrics
from surprise import (accuracy)
from sklearn.metrics import ndcg_score
# ! pip install recmetrics

# For FIM
import fim

# For Plotting
import matplotlib.pyplot as plt
import seaborn as sns
from itertools import cycle
import matplotlib.cm as cm
from plotly.subplots import make_subplots
import plotly.graph_objects as gob
from matplotlib import pyplot
import seaborn as sns

# For Deep Learning + SFIT
from statsmodels.stats.proportion import binom_test
from keras.initializers import glorot_normal
import keras.backend as K
from sklearn.model_selection import train_test_split
from sklearn.cluster import DBSCAN, SpectralClustering
from sklearn.cluster import AgglomerativeClustering, OPTICS, KMeans
from sklearn.metrics import accuracy_score, balanced_accuracy_score
from keras.models import Model, load_model
from keras.optimizers import Adam
from keras.layers import Input, Dense, Dropout, BatchNormalization
from keras.callbacks import EarlyStopping, ReduceLROnPlateau
from keras.utils import np_utils
import tensorflow as tf
from sfit import *

# Note: No Gap statistic to lessen runtime

def cluster_range(X, clusterer, k_start, k_stop, actual=None):
    '''
    Accepts the design matrix, the clustering object, the initial and final
    values to step through, and, optionally, actual labels and returns a
    dictionary of the cluster labels, cluster centers, internal validation
    values and, if actual labels is given, external validation values,
    for every k.
    '''
    ys = []
    centers = []
    inertias = []
    chs = []
    scs = []
    dbs = []
    gss = []
    gssds = []
    ps = []
    amis = []
    ars = []
    for k in range(k_start, k_stop+1):
        clusterer_k = clone(clusterer)
        clusterer_k.n_clusters = k
        clusterer_k.fit(X)
        y = clusterer_k.labels_
        ys.append(y)
        chs.append(calinski_harabasz_score(X, y))
        scs.append(silhouette_score(X, y, random_state=1337))
        dbs.append(davies_bouldin_score(X, y))
   
        centers.append(clusterer_k.cluster_centers_)
        inertias.append(clusterer_k.inertia_)
        
        if actual is not None:
            ps.append(purity(actual, y))
            amis.append(adjusted_mutual_info_score(actual, y))
            ars.append(adjusted_rand_score(actual, y))
            
    ret_dict = {'ys': ys, 'centers': centers, 'inertias': inertias,
                'chs': chs, 'scs': scs, 'dbs': dbs}
    
    if actual is not None:
        ret_dict['ps'] = ps
        ret_dict['amis'] = amis
        ret_dict['ars'] = ars
    
    return ret_dict

def plot_clusters(X, ys, centers):
    """Plot clusters given the design matrix and cluster labels"""
    k_max = len(ys) + 1
    k_mid = k_max//2 + 2
    fig, ax = plt.subplots(2, k_max//2, dpi=100, sharex=True, sharey=True, 
                           figsize=(16,9))
    for k, y in zip(range(2, k_max+1), ys):
        if k < k_mid:
            ax[0][k%k_mid-2].scatter(X[:,0], X[:,1], c=y, s=24, alpha=0.6,
                                     cmap='flare');
            ax[0][k%k_mid-2].set_title('$k=%d$ clusters'%k, size=18)
        else:
            ax[1][k%k_mid].scatter(X[:,0], X[:,1], c=y, s=24, alpha=0.6,
                                   cmap='flare');
            ax[1][k%k_mid].set_title('$k=%d$ clusters'%k, size=18)
    return ax

def plot_internal(inertias, chs, scs, dbs, gss=None, gssds=None):
    """Plot internal validation values"""
    fig, ax = plt.subplots()
    ks = np.arange(2, len(inertias)+2)
    ax.plot(ks, inertias, '-o', label='Sum of Squared Error',
            c='xkcd:chestnut')
    ax.plot(ks, chs, '-^', label='Calinski-Harabasz index',
            c='xkcd:peach')
    ax.set_xlabel('$k$ clusters', size=14)
    ax.set_ylabel('SSE/CH', size=14)
    lines, labels = ax.get_legend_handles_labels()
    ax2 = ax.twinx()
    ax2.plot(ks, scs, '-x', label='Silhouette coefficient',
             c='xkcd:dark pink')
    ax2.plot(ks, dbs, '-d', label='Davies-Bouldin index',
             c='xkcd:indigo')
    ax2.set_ylabel('Silhouette/DB', size=14)
    lines2, labels2 = ax2.get_legend_handles_labels()
    ax2.legend(lines+lines2, labels+labels2)
    return ax

def pooled_within_ssd(X, y, centroids, dist):
    """
    Compute pooled within-cluster sum of squares around the cluster mean
    """

    dist2 = []
    X = np.array(X)
    y = np.array(y)
    for x in np.unique(y):
        distance = [dist(i, centroids[y[index]])**2 for index,
                    i in enumerate(X) if y[index] == x]
        dist2.append(distance)
    
    w_k = sum([np.sum(i)/(2*len(i)) for i in dist2])
    return w_k

def cluster_range_kmedoids(X, k_start, k_stop, actual=None):
    """
    Accepts the design matrix, the initial and final values to step through,
    and, optionally, actual labels and returns a dictionary of the cluster
    labels, cluster centers, internal validation values, and external
    validation values if actual labels is given, for every k.
    """
    ys = []
    cs = []
    inertias = []
    chs = []
    scs = []
    dbs = []
    gss = []
    gssds = []
    ps = []
    amis = []
    ars = []
    for k in range(k_start, k_stop+1):
        clusterer_k = kmedoids(X, np.arange(k), ccore=True)
        clusterer_k.process()
        y = np.zeros(len(X), dtype=int)
        clusters = clusterer_k.get_clusters()
        for cluster, point in enumerate(clusters):
            y[point] = cluster
        centers = X[clusterer_k.get_medoids()]
        
        se = [euclidean(x, centers[y[i]]) ** 2 for i, x in enumerate(X)]
        sse = np.sum(se)
        ys.append(y)
        cs.append(centers)
        inertias.append(sse)
        chs.append(calinski_harabasz_score(X, y))
        scs.append(silhouette_score(X, y))
        dbs.append(davies_bouldin_score(X, y))

        if actual is not None:
                ps.append(purity(actual, y))
                amis.append(adjusted_mutual_info_score(actual, y))
                ars.append(adjusted_rand_score(actual, y))
                
    out = {
        'ys': ys,
        'centers': cs,
        'inertias': inertias,
        'chs': chs,
        'scs': scs,
        'dbs': dbs,
    }

    if actual is not None:
        out['ps'] = ps
        out['amis'] = amis
        out['ars'] = ars
    return out

def plot1(Z):
    '''
    Accepts a hierarchical clustering encoded as a linkage matrix and
    returns a truncated dendogram plot
    '''
    fig, ax = plt.subplots()
    dn = dendrogram(Z, ax=ax, truncate_mode='level', p=4)
    ax.set_ylabel(r'$\Delta$');
    
    return fig

def ave_nn_dist(n_neighbors, data):
    '''
    Accepts the n_neighbors to consider and the data then returns the average
    distance of a point to its  𝑘=1 nearest neighbor up to  𝑘=𝑛_𝑛𝑒𝑖𝑔ℎ𝑏𝑜𝑟𝑠 
    nearest neighbor as a sorted list.
    '''

    kNN = NearestNeighbors(n_neighbors=n_neighbors)
    kNN.fit(data)
    d = kNN.kneighbors()[0].mean(axis=1)
    
    return np.sort(d).tolist()

def plot_clusters_year(X, ys, centers):
    """Plot clusters given the design matrix and cluster labels"""
    k_max = len(ys) + 1
    k_mid = k_max//2 + 2
    fig, ax = plt.subplots(1, k_max-1, dpi=150, sharex=True, sharey=True, 
                           figsize=(7,4), subplot_kw=dict(aspect='equal'),
                           gridspec_kw=dict(wspace=0.01))
    for k, y in zip(range(2, k_max+1), ys):
            ax[k%k_mid-2].scatter(X[:,0], X[:,1], c=y, s=1, alpha=0.8);
            ax[k%k_mid-2].set_title('$k=%d$'%k)
    return ax

def plot_gmm(gmm, X, label=True, ax=None):
    ax = ax or plt.gca()
    labels = gmm.fit(X).predict(X)
    if label:
        ax.scatter(X[:, 0], X[:, 1], c=labels, s=40, cmap='viridis', zorder=2)
    else:
        ax.scatter(X[:, 0], X[:, 1], s=40, zorder=2)

    wfactor = 0.2 / gmm.weights.max()
    for pos, covar, w in zip(gmm.means, gmm.covariances, gmm.weights_):
        draw_ellipse(pos, covar, alpha=w * wfactor)
    plt.title("GMM with %d components"%len(gmm.means), fontsize=(20))
    plt.xlabel("U.A.")
    plt.ylabel("U.A.")


def SelBest(arr:list, X:int)->list:
    '''
    returns the set of X configurations with shorter distance
    '''
    dx=np.argsort(arr)[:X]
    return arr[dx]

def recsys(df_melt, rating_scale, method, user_list, df_name):
    '''
    Accepts a user rating dataframe, tuple of rating scale, 
    predition method, and list of user to give recommendations
    then returns a dataframe of the top 5 liked and recommended
    items per user based on the prediction method passed.
    '''
    df_rating = df_melt.pivot_table(values='rating', index='user_id',
                                    columns='recipe_id')
    reader = Reader(rating_scale=rating_scale)
    dataset = Dataset.load_from_df(df_melt, reader)
    method.fit(dataset.build_full_trainset())
    predictions = method.test(method.trainset.build_anti_testset())

    df_reco_user = pd.DataFrame()
    for user in user_list: # Get reco for N users
        top_n = []
        for uid, iid, true_r, est, _ in predictions:
            if uid == user:
                top_n.append((iid, est))
        
        # Get top N reco
        reco_list = sorted(top_n, key=lambda x: (-x[1], x[0]))[:10] 
        top = [i[0] for i in reco_list]
        reco = df_name.loc[top].reset_index(drop=True)\
                                .rename(columns={"name": user}).T
        x = df_name.loc[top].reset_index(drop=True)['name'].to_list()[:5]
        likes = df_rating.loc[user].dropna()\
                                    .sort_values(ascending=False).index
        y = df_name.loc[likes].reset_index(drop=True)['name'].to_list()[:5]
        assets = [user]
        Asset_feature = ['Likes', 'Recommends']
        index = ['Top1', 'Top2', 'Top3', 'Top4', 'Top5']

        midx = pd.MultiIndex.from_product([assets, Asset_feature])
        reco = pd.DataFrame([y,x], index=midx, columns=index).T
        df_reco_user = pd.concat([df_reco_user, reco], axis=1)

    return df_reco_user, predictions

def get_users_predictions(user_id, n, model):
    recommended_items = pd.DataFrame(model.loc[user_id])
    recommended_items.columns = ["predicted_rating"]
    recommended_items = recommended_items.sort_values('predicted_rating',
                                                      ascending=False)    
    recommended_items = recommended_items.head(n)
    return recommended_items.index.tolist()

def prepare_data(df, add_intercept=True):
    if add_intercept:
        df.insert(loc=0, column='intercept', value=1)
    msk = np.random.rand(len(df)) < 0.7
    train = df[msk]
    val_test = df[~msk]
    msk = np.random.rand(len(val_test)) < 0.66
    val = val_test[msk]
    test = val_test[~msk]
    print('Number of observations in train set: {0}'.format(train.shape[0]))
    print('Number of observations in validation set: {0}'\
          .format(val.shape[0]))
    print('Number of observations in test set: {0}'.format(test.shape[0]))
    print('Total number of observations: {0}'.format(train.shape[0] +
                                                     val.shape[0] +
                                                     test.shape[0]))
    X_train = train.iloc[:, 0:-1]
    Y_train = train.iloc[:, -1]
    X_val = val.iloc[:, 0:-1]
    Y_val = val.iloc[:, -1]
    X_test = test.iloc[:, 0:-1]
    Y_test = test.iloc[:, -1]
    return X_train, Y_train, X_val, Y_val, X_test, Y_test


def absolute_loss(predicted_y, true_y):
    """
    Compute the absolute value of the difference between two given
    real vectors of size N

    Parameters
    ----------
    predicted_y : numpy array of shape (N, 1)
        The predicted outcomes
    true_y : numpy array of shape (N, )
        The true outcomes

    Returns
    -------
    absolute_loss
        a numpy array of shape (N, )
    """
    return np.abs(true_y - np.squeeze(predicted_y))


def binary_cross_entropy_loss(predicted_y, true_y):
    """Compute the binary cross entropy loss between a vector of labels of
    size N and a vector of probabilities of same size

    Parameters
    ----------
    predicted_y : numpy array of shape (N, 1)
        The predicted probabilities
    true_y : numpy array of shape (N, )
        The true labels

    Returns
    -------
    binary_cross_entropy_loss
        a numpy array of shape (N, )
    """
    x = np.squeeze(predicted_y)
    return -np.log(x)*true_y - np.log(1 - x)*(1 - true_y)


def categorical_cross_entropy_loss(predicted_y, true_y):
    """Compute the categorical cross entropy loss between a vector of
    labels of size (N,C) and a vector of probabilities of same size

    Parameters
    ----------
    predicted_y : numpy array of shape (N, 1)
        The predicted probabilities
    true_y : numpy array of shape (N, )
        The true labels

    Returns
    -------
    categorical_cross_entropy_loss
        a numpy array of shape (N, )
    """
    return -np.sum(np.multiply(np.log(np.maximum(predicted_y, 1e-6)),
                               true_y), axis=1)


def compute_errors_first(loss, model, n, p, x, y):
    x_intercept = np.copy(x)
    x_intercept[:, range(1, p)] = 0
    predicted_y = model.predict(x_intercept)
    errors = np.zeros((n, p))
    errors[:, 0] = loss(predicted_y, y)
    for j in range(1, p):
        x_j = np.copy(x)
        indices = [i for i in range(1, p) if i != j]
        x_j[:, indices] = 0
        predicted_y_j = model.predict(x_j)
        errors[:, j] = loss(predicted_y_j, y)
    return errors


def sign_test_first(alpha, beta, errors, n, p):
    c_1 = {}
    s_1 = []
    u_1 = []
    p_values = np.zeros(p - 1)
    for j in range(1, p):
        delta_j = (1 - beta) * errors[:, 0] - errors[:, j]
        n_j = np.sum(delta_j > 0)
        p_values[j - 1] = binom_test(n_j, n, 0.5, 'larger')
        if binom_test(n_j, n, 0.5, 'larger') < alpha:
            q = norm.ppf(1 - alpha / 2, 0, 1)
            lower = int(np.floor((n + 1) / 2 - q * np.sqrt(n) / 2))
            upper = int(np.ceil((n + 1) / 2 + q * np.sqrt(n) / 2))
            ordered_delta_j = np.sort(delta_j)
            s_1.append(j)
            c_1[j] = (np.median(delta_j), (ordered_delta_j[lower],
                                           ordered_delta_j[upper]))
        else:
            u_1.append(j)
    return c_1, p_values, s_1, u_1


def sfit_first_order(model, loss, x, y, alpha, beta=None, verbose=0):
    """
    Compute the first-order SFIT method to test what are the
    first-order significant variables within x toward the prediction
    of y as learned by the model.

    Parameters
    ----------
    model :
        A predictive model that can return predicted y given x by calling its
        method model.predict
    loss : function
        Function that computes the pointwise loss between 2 numpy arrays of
        outcomes, its first argument should be the predicted outcomes and its
        second should be the true ones.
    x: numpy array of shape (N, p)
        Input data used to perform the tests
    y: numpy array of shape (N, )
        True outcomes
    alpha : float
        Significance level of the test
    beta: float
        Regularization amount of the test. If set to None, the optimal beta
        parameter is found using the randomization procedure described in the
        original SFIT paper (works only when model is a Keras model).
    verbose: boolean
        The summary of the test procedure is printed if true (default)
        but no printing if false.
    Returns
    -------
    s_1 : list
        the list of first-order significant variables (indexed from 1 to p)
    c_1 : dictionary
        dictionary whose keys are the first-order significant variables ;
        for each key, its value is a tuple whose first element is the test
        statistic value and second element is its (1 - alpha)%
        confidence interval
    u_1 : list
        the list of first-order non-significant variables
        (indexed from 1 to p)
    p_values : numpy array of shape (p, )
        array containing the p-values associated with each variables.
    second_order_significance : boolean
        If true, indicates the presence of significant second-order effects
        which suggests to use second-order SFIT.
    opt_beta : float
        Equals to beta if beta is passed as argument, otherwise equals
        to the optimal beta found from randomization
        procedure.
    """
    (n, p) = x.shape
    assert np.array_equal(x[:, 0], np.ones(n)), 'As an intercept, all the '\
                                            'entries of the first column '\
                                            'of x should be equal to 1.'
    errors = compute_errors_first(loss, model, n, p, x, y)
    if beta is not None:
        opt_beta = beta
        c_1, p_values, s_1, u_1 = sign_test_first(alpha, beta, errors, n, p)
    else:
        print('Compute optimal beta')
        beta = 1e-7
        mean_nr_significants = 1
        initial_weights = model.get_weights()
        k_eval = lambda placeholder: placeholder.eval(session=K.get_session())
        nr_simulations = 20
        x_val, x_test, y_val, y_test = train_test_split(x, y, test_size=0.5)
        n_val = x_val.shape[0]
        n_test = x_test.shape[0]
        while mean_nr_significants > alpha:
            beta = beta*10
            nr_significants_per_sim = np.zeros(nr_simulations)
            for i in range(nr_simulations):
                new_weights = [k_eval(glorot_normal()(w.shape)) \
                               for w in initial_weights]
                model.set_weights(new_weights)
                errors_i = compute_errors_first(loss, model, n_val, p,
                                                x_val, y_val)
                c_1_i, p_values_i, s_1_i, u_1_i = sign_test_first(alpha,
                                                                  beta,
                                                                  errors_i,
                                                                  n_val, p)
                nr_significants_per_sim[i] = len(s_1_i)/p
            mean_nr_significants = np.mean(nr_significants_per_sim)
        opt_beta = beta
        print('Optimal beta found: {0}'.format(opt_beta))
        model.set_weights(initial_weights)
        errors = compute_errors_first(loss, model, n_test, p, x_test, y_test)
        c_1, p_values, s_1, u_1 = sign_test_first(alpha, beta, errors,
                                                  n_test, p)
    if verbose:
        print('Summary of first-order SFIT\n'
              '------------------------------------------------\n'
              'First-order significant variables:')
        for key in c_1.keys():
            print('- Variable {0}:'.format(key))
            print('\t Median: {0}'.format(np.round(c_1[key][0], 3)))
            print('\t {0}% confidence interval: {1}'\
                  .format(int(100*(1-alpha)), np.round(c_1[key][1], 3)))
        print('------------------------------------------------\n'
              'First-order non-significant variables: {0}'.format(u_1))

    # Test for presence of any second-order significance:
    second_order_significance = False
    x_first = np.copy(x)
    indices = [i for i in range(1, p) if i in u_1]
    x_first[:, indices] = 0
    predicted_y_first = model.predict(x_first)
    model_first_errors = loss(predicted_y_first, y)
    predicted_y_all = model.predict(x)
    model_all_errors = loss(predicted_y_all, y)
    delta_j = model_first_errors - model_all_errors
    n_j = np.sum(delta_j > 0)
    if binom_test(n_j, n, 0.5, 'larger') < alpha:
        second_order_significance = True
        if verbose:
            print('------------------------------------------------\n'
                  'There are some significant second-order variables: '\
                  'recommended to run second-order SFIT.')
    else:
        if verbose:
            print('------------------------------------------------\n'
                  'There are no significant second-order variables.')
    if verbose:
        print('------------------------------------------------ \n'
              '------------------------------------------------ \n')
    return s_1, c_1, u_1, p_values, second_order_significance, opt_beta


def sfit_second_order(model, loss, alpha, beta, x, y, s_1, u_1, verbose=0):
    """Compute the second-order SFIT method to test what are the second-order
    significant variables within x toward the prediction of y as learned
    by the model.

    Parameters
    ----------
    model :
        A predictive model that can return predicted y given x by calling
        its method model.predict
    loss : function
        Function that computes the pointwise loss between 2 numpy arrays of
        outcomes, its first argument should be the predicted outcomes and
        its second should be the true ones.
    alpha : float
        Significance level of the test
    beta: float
        Regularization amount of the test
    x: numpy array of shape (N, p)
        Input data used to perform the tests
    y: numpy array of shape (N, )
        True outcomes
    s_1 : list
        the list of first-order significant variables as returned by
        sfit_first_order
    u_1 : list
        the list of first-order non-significant variables as returned by
        sfit_first_order
    verbose: boolean
        The summary of the test procedure is printed if true (default) but
        no printing if false.

    Returns
    -------
    s_2: list
        the list of second-order significant variables (indexed from 1 to p)
    c_2 : dictionary
        dictionary whose keys are the second-order significant pairs pf
        variables ; for each key, its value is a tuple whose first element
        is the test statistic value and second element is its (1 - alpha)%
        confidence interval
    u_2 : list
        the list of second-order non-significant variables
        (indexed from 1 to p)
    third_order_significance : boolean
        If true, indicates the presence of significant third-order effects.
    """
    s_2 = []
    c_2 = {}
    u_2 = []
    (n, p) = x.shape
    assert np.array_equal(x[:, 0], np.ones(n)), 'As an intercept, all the '\
    'entries of the first column of x should be equal to 1.'
    x_intercept = np.copy(x)
    x_intercept[:, range(1, p)] = 0
    predicted_y_intercept = model.predict(x_intercept)
    baseline_model_errors = loss(predicted_y_intercept, y)
    for j in u_1:
        for k in range(1, p):
            x_jk = np.copy(x)
            indices = [i for i in range(1, p) if (i != j and i != k)]
            x_jk[:, indices] = 0
            predicted_y_jk = model.predict(x_jk)
            model_jk_errors = loss(predicted_y_jk, y)
            if k in s_1:
                x_k = np.copy(x)
                indices = [i for i in range(1, p) if i != k]
                x_k[:, indices] = 0
                predicted_y_k = model.predict(x_k)
                model_k_errors = loss(predicted_y_k, y)
                delta_jk = (1-beta)*model_k_errors - model_jk_errors
            else:
                delta_jk = (1-beta)*baseline_model_errors - model_jk_errors
            n_jk = np.sum(delta_jk > 0)
            if binom_test(n_jk, n, 0.5, 'larger') < alpha:
                q = norm.ppf(1-alpha/2, 0, 1)
                lower = int(np.floor((n+1)/2 - q*np.sqrt(n)/2))
                upper = int(np.ceil((n+1)/2 + q*np.sqrt(n)/2))
                ordered_delta_jk = np.sort(delta_jk)
                if j not in s_2:
                    s_2.append(j)
                c_2[(j, k)] = (np.median(delta_jk), (ordered_delta_jk[lower],
                                                     ordered_delta_jk[upper]))
        if j not in s_2:
            u_2.append(j)
    if verbose:
        print('Summary of second-order SFIT\n'
              '------------------------------------------------\n'
              'Second-order significant variables:')
        for key in c_2.keys():
            print('- Variable {0}:'.format(key))
            print('\t Median: {0}'.format(np.round(c_2[key][0], 3)))
            print('\t {0}% confidence interval: {1}'\
                  .format(int(100*(1-alpha)), np.round(c_2[key][1], 3)))
        print('------------------------------------------------\n'
              'Second-order non-significant variables: {0}'.format(u_2))

    # Test for presence of any third-order significance:
    x_second = np.copy(x)
    indices = [i for i in range(1, p) if i in u_2]
    x_second[:, indices] = 0
    predicted_y_second = model.predict(x_second)
    model_second_errors = loss(predicted_y_second, y)
    predicted_y_all = model.predict(x)
    model_all_errors = loss(predicted_y_all, y)
    delta_j = model_second_errors - model_all_errors
    n_j = np.sum(delta_j > 0)
    third_order_significance = False
    if binom_test(n_j, n, 0.5, 'larger') < 0.05:
        third_order_significance = True
        if verbose:
            print('------------------------------------------------\n'
                  'There are some significant third-order variables: '\
                  'recommended to run third-order SFIT.')
    else:
        if verbose:
            print('------------------------------------------------\n'
                  'There are no significant third-order variables.')
    if verbose:
        print('------------------------------------------------ \n'
              '------------------------------------------------ \n')

    return s_2, c_2, u_2, third_order_significance

pd.read_csv('./RAW_recipes.csv',
            index_col='id').head(5)

pd.read_csv('./RAW_interactions.csv', index_col='user_id').head(5)

# filter foods tagged as dessert
df = pd.read_csv('./RAW_recipes.csv', index_col='id')
cols = ['minutes', 'tags', 'nutrition', 'n_steps', 'ingredients']
df_dessert = df.loc[df['tags'].str.contains("dessert", case=False)][cols]
df_dessert = df_dessert.loc[(df['minutes'] <= 90) & (df['n_steps'] <= 7)]

# Store id vs name of dessert
df_name = df.loc[df['tags'].str.contains("dessert", case=False)][['name']]

# Extracting columns from nutrition list
df_dessert['nutrition'] = df_dessert['nutrition'].apply(ast.literal_eval) 
df_dessert['ingredients'] = df_dessert['ingredients'].apply(ast.literal_eval)
df_dessert['ingredients'] = [', '.join(map(str, l)) \
                             for l in df_dessert['ingredients']]
col_nutri = ['calories', 'total_fat', 'sugar_pvd', 'sodium_pvd',
             'protein_pvd', 'sat_fat_pvd', 'carbs_pvd']
df_nutri = pd.DataFrame(df_dessert['nutrition'].to_list(), columns=col_nutri,
                        index=df_dessert.index)
df_dessert2 = df_dessert.join(df_nutri).drop(['nutrition', 'tags'], axis=1)

# Extracting columns from ingredient's list then onehot encode ingredients
tfidf = TfidfVectorizer(
    analyzer=lambda d: d.split(', '),
    stop_words='english',
    max_df=0.9,
    min_df=0.03
)
ingredients = tfidf.fit_transform(df_dessert2['ingredients'])
df_recipe = pd.DataFrame(ingredients.todense(),
                         columns=tfidf.get_feature_names_out(),
                         index=df_dessert.index)
df_dessert3 = df_dessert2.join(df_recipe).drop('ingredients', axis=1)
df_dessert_all = df_dessert3.drop(index=46942)

# Store recipe id's from our items profile
food_list = df_dessert_all.index.tolist()

# Load Raw Interactions Dataset
file_path = './RAW_interactions.csv'
df_pp_recipes = pd.read_csv(file_path)

# Parse dates from string to datetime
df_pp_recipes['date'] = pd.to_datetime(df_pp_recipes['date'],
                                       format='%Y-%m-%d')

# Filter users who rated from 2012 to 2018
filtered_df = df_pp_recipes.loc[(df_pp_recipes['date'] >= '2008-01-01')]
# filtered_df = df_pp_recipes.copy()

# Filter users who rated for items in the items profile
filtered_df = filtered_df[filtered_df['recipe_id'].isin(food_list)]
filtered_df = filtered_df[filtered_df.groupby('user_id')['user_id']\
                          .transform('size') >= 5]

# Create utility matrix
df_user_ratings = pd.pivot_table(filtered_df,
                                values='rating',
                                index='user_id',
                                columns='recipe_id')

df_user_ratings

# Filtering recipes matrix to retain only items rated by users
rated_items = (df_user_ratings.columns.tolist())
df_dessert_final = df_dessert_all.loc[df_dessert_all.index.isin(rated_items)]
df_dessert_final

text_list = df_dessert['ingredients'].to_list()
text = ''.join(text_list)

mask = np.array(Image.open("chef.png"))

wc = WordCloud(width = 3000, height = 2000, random_state=1,
               max_words = 100, max_font_size = 500,
               background_color='white', colormap='flare',
               stopwords = STOPWORDS, mask=mask)
wc.generate(text)

plt.figure( figsize=(16,8))
plt.imshow(wc, interpolation="None")
plt.axis('off')
plt.show()

text_list = df_name['name'].to_list()
text = ''.join(text_list)

mask = np.array(Image.open("chef.png"))
STOPWORDS = ['S', 'and', 'or', 'with']
wc = WordCloud(width = 3000, height = 2000, random_state=1,
               max_words = 100, max_font_size = 500,
               background_color='white', colormap='flare',
               stopwords = STOPWORDS, mask=mask)
wc.generate(text)

plt.figure( figsize=(16,8))
plt.imshow(wc, interpolation="None")
plt.axis('off')
plt.show()

a4_dims = (12, 6)
fig, ax = pyplot.subplots(figsize=a4_dims)
ax = sns.histplot(data=filtered_df,
             x='rating',
             kde=True,
             color='navy',)
 
plt.xlabel("Rating", fontsize = 20)
plt.ylabel("Count", fontsize = 20)
plt.tick_params(axis='both', which='major', labelsize=14)
plt.title("Histogram of Dessert Rating", fontsize = 24);

count = filtered_df.groupby('recipe_id')['rating'].count()\
                                                .sort_values(ascending=False)

a4_dims = (12, 6)
fig, ax = pyplot.subplots(figsize=a4_dims)
ax = sns.histplot(data=count.reset_index(),
             x='rating',
             kde=True,
             color='navy',)
 
plt.xlabel("Rating Count", fontsize = 20)
plt.ylabel("Count", fontsize = 20)
plt.tick_params(axis='both', which='major', labelsize=14)
plt.title("Histogram of Dessert Rating Count", fontsize = 24);

labels = ["Kittencal's Choco Frosting",
          'Vanilla  Buttercream Frosting',
          'My No Roll Pie Crust',
          'Reeses Squares',
          "Edna's Apple Crumble"]
colors = ['chocolate', 'olive', 'tan', 'orange', 'red']
values = count[:5].values

height = 0.9
plt.figure(figsize=(12, 6))
plt.barh(y=labels, width=values, height=height, color=colors, align='center')
plt.xlim(0, max(values) * 1.2)
plt.ylim(-0.5, len(labels) - 0.5)
for i, v in enumerate(values):
    plt.text(v+5, i, str(v), color=colors[i], fontweight='bold')
plt.title('Highest Number of Ratings of Desserts',
          fontdict={'fontsize': 24,
                  'verticalalignment': 'baseline'})
plt.xlabel("Number of Rating", fontsize=16)
plt.ylabel("Dessert Name", fontsize=16)
plt.tick_params(axis='both', which='major', labelsize=14)
plt.show()

# 
sns.set_style({'axes.grid' : False})
palette = sns.color_palette('flare')
sns.set_style('darkgrid')
# sns.set(rc={'figure.figsize': (20, 6)})

# Perform SVD (untruncated)
df_user_rating_0 = df_user_ratings.fillna(-1)
df_user_rating_0

# Scaling Dataset

standard_scaler = StandardScaler()
X_user = standard_scaler.fit_transform(df_user_rating_0.T).T

n_components = len(df_user_rating_0.columns)-1
svd_user = TruncatedSVD(n_components=n_components, random_state=1337)
X_svd_user = svd_user.fit_transform(X_user)
variance_explained_user = svd_user.explained_variance_ratio_
p_svd_user = svd_user.components_

# Plotting the cumulative variance explained
fig, ax = plt.subplots(figsize=(12, 6))
plt.plot(range(1, len(variance_explained_user)+1), 
         variance_explained_user.cumsum(), '-', c='xkcd:chestnut',
         linewidth=2.5)
plt.ylim(0,1)
plt.xlabel('number of singular vectors')
plt.ylabel('cumulative variance explained')
plt.axhline(y=0.8, color='orange', linestyle='--', linewidth=2)
plt.axvline(x=np.where(variance_explained_user.cumsum() > 0.81)[0][0],
            color='orange', linestyle='--', linewidth=2);

# Perform TruncatedSVD

component_user = np.where(variance_explained_user.cumsum() > 0.8)[0][0]
svd_user = TruncatedSVD(n_components=component_user, random_state=1337)
X_svd_user = svd_user.fit_transform(X_user)

# Plotting the first two Singular Vectors
fig, ax = plt.subplots(figsize=(12, 6))
plt.scatter(X_svd_user[:,0], X_svd_user[:,1], c='steelblue', alpha=0.6)
plt.xlabel('SV1', size=20)
plt.ylabel('SV2', size=20);

feature_names = df_user_ratings.columns
for i in range(5):
    fig, ax = plt.subplots(figsize=(12, 6))
    order = np.argsort(np.abs(p_svd_user[:, i]))[-20:]
    ax.barh([df_name.loc[feature_names[o]][0] for o in order],
            p_svd_user[order, i], color='steelblue')
    ax.set_title(f'SV{i+1}', size=20)

# Perform SVD (untruncated)

# Scaling Dataset
standard_scaler = StandardScaler()
X_item = standard_scaler.fit_transform(df_dessert_final)

n_components = len(df_dessert_final.columns)-1
svd_item = TruncatedSVD(n_components=n_components, random_state=1337)
X_svd_item = svd_item.fit_transform(X_item)
variance_explained_item = svd_item.explained_variance_ratio_
p_svd_item = svd_item.components_

# Plotting the cumulative variance explained
fig, ax = plt.subplots(figsize=(12, 6))
plt.plot(range(1, len(variance_explained_item)+1), 
         variance_explained_item.cumsum(), '-', c='xkcd:chestnut',
         linewidth=2.5)
plt.ylim(0,1)
plt.xlabel('number of singular vectors')
plt.ylabel('cumulative variance explained')
plt.axhline(y = 0.8, color = 'orange', linestyle = '--')
plt.axvline(x = np.where(variance_explained_item.cumsum() > 0.81)[0][0],
            color = 'orange', linestyle = '--');

# Perform TruncatedSVD
component_item = np.where(variance_explained_item.cumsum() > 0.8)[0][0]
svd_item = TruncatedSVD(n_components=component_item, random_state=1337)
X_svd_item = svd_item.fit_transform(X_item)

# Plotting the first two Singular Vectors
fig, ax = plt.subplots(figsize=(12, 6))
plt.scatter(X_svd_item[:,0], X_svd_item[:,1], c='xkcd:chocolate brown',
            alpha=0.6)
plt.xlabel('singular vector 1')
plt.ylabel('singular vector 2');